{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "import pymc3 as pm\n",
    "import theano.tensor as tt\n",
    "import scipy.stats as stats\n",
    "from sklearn.metrics import confusion_matrix,roc_auc_score,roc_curve\n",
    "from mcmc import mcmc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 数据生成与参数估计\n",
    "# P值检验和指标获取\n",
    "# 指标保存\n",
    "# 多线程"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=500\n",
    "I=10\n",
    "ab_p = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 数据生成\n",
    "lambda_i = np.random.normal(4,0.25,size=(1,I))\n",
    "phi_i = np.random.normal(1,0.17*0.5,size=(1,I))\n",
    "zeta_p = np.random.normal(0,1,size=(N,1))\n",
    "zeta_p_speedness = np.random.normal(1.5,1,size=(N,1))\n",
    "epsilon = 1/np.random.uniform(1,5,size=(1,I))\n",
    "rt = np.random.lognormal(phi_i*(lambda_i-zeta_p),epsilon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 数据模拟\n",
    "## random 随机作答\n",
    "random_flag = np.random.binomial(1,ab_p,(N,1))\n",
    "random_rt = np.random.normal(4,1,size=(N,I))*random_flag+rt*(1-random_flag)\n",
    "## speedness 加速作答\n",
    "speedness_flag = np.random.binomial(1,ab_p,(N,1))*np.random.binomial(1,0.5,(N,I))\n",
    "speedness_rt = np.random.lognormal(phi_i*(lambda_i-zeta_p_speedness),epsilon**0.5)*speedness_flag+\\\n",
    "                rt*(1-speedness_flag)\n",
    "## extrme 极端作答\n",
    "extrme_flag = np.random.binomial(1,ab_p,(N,1))*np.random.randint(0,I,(N,1))==np.arange(I)\n",
    "extrme_rt = lambda_i.max()*2*extrme_flag+rt*(1-extrme_flag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def l_s(z_square,I):\n",
    "    c1 = (z_square.sum(axis=1)/I)**(1/3)\n",
    "    c2 = (1-2/(9*I))\n",
    "    c3 = (2/(9*I))**0.5\n",
    "    return (c1-c2)/c3\n",
    "def l_0(z_square,epsilon):\n",
    "    c1=z_square\n",
    "    c2 = tt.log(2*np.pi*(epsilon**2))\n",
    "    return (c1+c2).sum(axis=1)\n",
    "def l_z(z_square,I):\n",
    "    c1 = z_square.sum(axis=1)\n",
    "    c2 = I\n",
    "    c3 = (2*I)**0.5\n",
    "    return (c1-c2)/c3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_mat(l_0,l_z,l_s,I):\n",
    "    z_p = stats.norm(0,1).cdf(l_z).mean(axis=0).mean(axis=0)\n",
    "    o_p = stats.chi2(I-1).cdf(l_0).mean(axis=0).mean(axis=0)\n",
    "    s_p = stats.norm(0,1).cdf(l_s).mean(axis=0).mean(axis=0)\n",
    "    return o_p,z_p,s_p\n",
    "def get_specific_sensitive(confusion_mat):\n",
    "    indicator = confusion_mat*np.eye(2)/confusion_mat.sum(axis=1)\n",
    "    # 敏感性 特异性\n",
    "    return indicator[0,0],indicator[1,1]\n",
    "def p_to_flag_pre(p_values,alpha):\n",
    "    return [p>(1-alpha) for p in p_values]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 4 jobs)\n",
      "NUTS: [epsilon, phi, phi_scale, phi_mean, lambda, lambda_scale, lambda_mean, zeta]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='20000' class='' max='20000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [20000/20000 02:47<00:00 Sampling 2 chains, 146 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 5_000 tune and 5_000 draw iterations (10_000 + 10_000 draws total) took 168 seconds.\n",
      "There were 146 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "The acceptance probability does not match the target. It is 0.6387110219908021, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "The estimated number of effective samples is smaller than 200 for some parameters.\n"
     ]
    }
   ],
   "source": [
    "\n",
    "def l_s(z_square,I):\n",
    "    c1 = (z_square.sum(axis=1)/I)**(1/3)\n",
    "    c2 = (1-2/(9*I))\n",
    "    c3 = (2/(9*I))**0.5\n",
    "    return (c1-c2)/c3\n",
    "def l_0(z_square,epsilon):\n",
    "    c1=z_square\n",
    "    c2 = tt.log(2*np.pi*(epsilon**2))\n",
    "    return (c1+c2).sum(axis=1)\n",
    "def l_z(z_square,I):\n",
    "    c1 = z_square.sum(axis=1)\n",
    "    c2 = I\n",
    "    c3 = (2*I)**0.5\n",
    "    return (c1-c2)/c3\n",
    "with pm.Model() as model:\n",
    "    zeta_pre = pm.Normal(\"zeta\",0,1,shape=(N,1))\n",
    "    lambda_mean = pm.Normal(\"lambda_mean\",4,2)\n",
    "    lambda_scale = pm.InverseGamma(\"lambda_scale\",1,1)\n",
    "    lambda_pre = pm.Normal(\"lambda\",lambda_mean,lambda_scale,shape=(1,I))\n",
    "    phi_mean = pm.Normal(\"phi_mean\",1,2)\n",
    "    phi_scale = pm.InverseGamma(\"phi_scale\",1,1)\n",
    "    phi_pre = pm.Normal(\"phi\",phi_mean,phi_scale,shape=(1,I))\n",
    "    epsilon_pre = pm.InverseGamma(\"epsilon\",1,1,shape=(1,I))\n",
    "    mu_ip = phi_pre*(lambda_pre-zeta_pre)\n",
    "    z_squre = ((tt.log(rt)-mu_ip)/epsilon_pre)**2 #  ((np.log(rt)-phi_i*(lambda_i-zeta_p))/epsilon)**2\n",
    "    l_0 = l_0(z_squre,epsilon_pre)\n",
    "    l_z = l_z(z_squre,I)\n",
    "    l_s = l_s(z_squre,I)\n",
    "    pm.Deterministic(\"z_squre\",z_squre)\n",
    "    pm.Deterministic(\"l_0\",l_0)\n",
    "    pm.Deterministic(\"l_z\",l_z)\n",
    "    pm.Deterministic(\"l_s\",l_s)\n",
    "    pm.Lognormal(\"T\",mu_ip,epsilon_pre,observed=random_rt)\n",
    "    trace = pm.sample(5000,tune=5000,chains=2,\n",
    "                        return_inferencedata=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "o_p,z_p,s_p = p_mat(trace.posterior[\"l_0\"],trace.posterior[\"l_z\"],trace.posterior[\"l_s\"],I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas = [0.2,0.05,0.01]\n",
    "indicatiors=[ get_specific_sensitive(confusion_matrix(random_flag,p_to_flag_pre(o_p,alpha),normalize=\"true\",labels=[1,0])) for alpha in alphas]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#3.保存\n",
    "def save(data):\n",
    "    import pickle\n",
    "    import time \n",
    "    while True:\n",
    "        try:\n",
    "            with open(\"./indicator.pkl\",\"rb\") as f:\n",
    "                indicator_saved = pickle.load(f)\n",
    "            indicator_saved.append(data)\n",
    "            with open(\"./indicator.pkl\",\"rb\") as f:\n",
    "                pickle.dump(indicator_saved,f) \n",
    "            break\n",
    "        except Exception as e:\n",
    "            print(e)\n",
    "            time.sleep(1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.多线程\n",
    "import threading\n",
    "class ModelRun(threading.Thread):  \n",
    "    def __init__(self, func, args=()):  \n",
    "        super(ModelRun, self).__init__()  \n",
    "        self.func = func  \n",
    "        self.args = args  \n",
    "  \n",
    "    def run(self):  \n",
    "        self.result = self.func(*self.args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([478,  22])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "confusion_matrix(random_flag,s_p>0.95).sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 5000, 500)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.norm(0,1).cdf(trace.posterior[\"l_z\"]).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAARfklEQVR4nO3dfaxkdX3H8fenrJqCtUD2QhHQxRYfgPhAb6nWxhApgYi62GqLjbpB7MZGrW2a1qVNSxNjs42Nsa1Vs1HKmloo8YmtT5WstaRNBBdBBBaECoWVLXstrbbaqIvf/jEHGS73cufO05372/cr2cycM2fmfPfcuZ/5zu883FQVkqS2/NhaFyBJGj/DXZIaZLhLUoMMd0lqkOEuSQ3asNYFAGzcuLE2bdq01mVI0rpy/fXXf7Oq5pZ6bCbCfdOmTezZs2ety5CkdSXJvy/32IrDMkkuTXIgyc19896Z5LYkNyX5eJIj+x67OMmdSW5Pcs7I1UuSVm2QMffLgHMXzbsaOK2qng18DbgYIMkpwAXAqd1z3pvksLFVK0kayIrhXlXXAA8smve5qjrYTX4ROKG7vxm4oqq+V1V3AXcCZ4yxXknSAMZxtMzrgc90948H7u17bF8371GSbE2yJ8mehYWFMZQhSXrISOGe5A+Bg8CHH5q1xGJLXrymqnZU1XxVzc/NLbmzV5I0pKGPlkmyBXgpcFY9fPWxfcCJfYudANw3fHmSpGEM1bknORd4G/Dyqvpu30O7gAuSPCHJScDJwHWjlylJWo0VO/cklwNnAhuT7AMuoXd0zBOAq5MAfLGq3lhVtyS5EriV3nDNm6rqwUkVL0laWmbheu7z8/PlSUyStDpJrq+q+aUe89oyh7hN2z611iVImgDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGe8O84qN06DLcJalBhnuDNm371KO6drt46dBiuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMM98Z4JqokGCDck1ya5ECSm/vmHZ3k6iR3dLdH9T12cZI7k9ye5JxJFS5JWt4gnftlwLmL5m0DdlfVycDubpokpwAXAKd2z3lvksPGVq0ewS5d0nJWDPequgZ4YNHszcDO7v5O4Py++VdU1feq6i7gTuCM8ZQqSRrUsGPux1bVfoDu9phu/vHAvX3L7evmPUqSrUn2JNmzsLAwZBmHhqWu8DhM1z7s8yStP+PeoZol5tVSC1bVjqqar6r5ubm5MZchSYe2YcP9/iTHAXS3B7r5+4AT+5Y7Abhv+PK0ErtxSUvZMOTzdgFbgO3d7VV98/8uybuAJwMnA9eNWqRWZsBL6rdiuCe5HDgT2JhkH3AJvVC/MslFwD3AqwCq6pYkVwK3AgeBN1XVgxOqXZK0jBXDvapevcxDZy2z/DuAd4xSlCRpNJ6hKkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkho07ElMWiPjOlnJk56kttm5S1KDDPcZM41rxdi1S+0z3CWpQY65N261XfpDy9+9/bxJlCNpSuzcJalBhrskNchwl6QGGe6S1CDDXT/iIZJSOwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXUvysEhpfTPcJalBXhVSdulSg+zcJalBI3XuSX4HeANQwFeBC4HDgb8HNgF3A79aVf81UpWyu5a0KkN37kmOB34LmK+q04DDgAuAbcDuqjoZ2N1NS5KmaNRhmQ3AjyfZQK9jvw/YDOzsHt8JnD/iOiRJqzR0uFfVN4A/B+4B9gPfqqrPAcdW1f5umf3AMUs9P8nWJHuS7FlYWBi2DK2BafwRb0mjGWVY5ih6XfpJwJOBI5K8ZtDnV9WOqpqvqvm5ublhy5AkLWGUYZlfAu6qqoWq+gHwMeAXgPuTHAfQ3R4YvUxJ0mqMEu73AM9PcniSAGcBe4FdwJZumS3AVaOVKElaraEPhayqa5N8BPgycBC4AdgBPBG4MslF9D4AXjWOQiVJgxvpOPequgS4ZNHs79Hr4iVJa8QzVCWpQYa7JDXIcJekBnlVyBnmiUKShmXnLkkNMtw1EC85IK0vhrskNchwnwF2xZLGzXDXsvzQkdYvw12SGmS4z5BZ7ZJntS5JyzPcJalBhvsaW64rtluWNArDXZIaZLhLUoMMd0lqkOEuSQ0y3DU0d/pKs8twl6QGGe6S1CDDXZIaZLhrVRxnl9YHw12SGmS4S1KDDHdJapDhLkkNGinckxyZ5CNJbkuyN8kLkhyd5Ookd3S3R42rWM0e/1qTNJtG7dz/AvhsVT0TeA6wF9gG7K6qk4Hd3bQkaYqGDvckTwJeBHwQoKq+X1X/DWwGdnaL7QTOH61ESdJqjdK5Pw1YAP4myQ1JPpDkCODYqtoP0N0eM4Y6JUmrMEq4bwBOB95XVc8DvsMqhmCSbE2yJ8mehYWFEcqQJC02SrjvA/ZV1bXd9Efohf39SY4D6G4PLPXkqtpRVfNVNT83NzdCGZKkxYYO96r6D+DeJM/oZp0F3ArsArZ087YAV41UoSRp1UY9WuYtwIeT3AQ8F/hTYDtwdpI7gLO7afVp9dBBD4uUZseGUZ5cVTcC80s8dNYorytJGo1nqEpSgwx3SWqQ4S5JDTLcJalBhrskNchw11is9hBID5uUJstwl6QGjXScu4Zn1yppkuzcJalBhrskNchwl6QGOeY+RY6zS5oWO3dJapDhLkkNMtw1dp6gJK09w12SGmS4S1KDDHdJapDhLkkNMtw1Me5UldaO4S5JDTLcNVEeFimtDcNdkhpkuE+J3aukaTLcJalBXhVywuzYJa0FO3dNjTtXpekZOdyTHJbkhiSf7KaPTnJ1kju626NGL1OStBrj6NzfCuztm94G7K6qk4Hd3bQkaYpGCvckJwDnAR/om70Z2Nnd3wmcP8o6JEmrN2rn/m7g94Ef9s07tqr2A3S3xyz1xCRbk+xJsmdhYWHEMiRJ/YYO9yQvBQ5U1fXDPL+qdlTVfFXNz83NDVuGJGkJo3TuLwRenuRu4ArgxUn+Frg/yXEA3e2BkatUU5Y7YsYjaaTxGTrcq+riqjqhqjYBFwCfr6rXALuALd1iW4CrRq5SkrQqkziJaTtwZZKLgHuAV01gHWqIHbs0fmMJ96r6AvCF7v5/AmeN43UlScPxDFVJapDhLkkNMtwlqUGG+wR4gSxJa81wl6QGGe6aitWeuOQ3H2k0hrvWheWGuhwCk5ZmuEtSgwz3CbKjHD+3qTQYw12SGmS4S1KDDHdJatAkrgopDc0xdWk8DHfNLINeGp7DMpLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLc1QwvIiY9zHCXpAYZ7pLUIMNd6844h18cxlGrDHdJapDhPibuzJsd/T8Hfy46VBnuktSgocM9yYlJ/inJ3iS3JHlrN//oJFcnuaO7PWp85UoPG3dHboevlozSuR8EfreqngU8H3hTklOAbcDuqjoZ2N1NS5KmaOhwr6r9VfXl7v7/AHuB44HNwM5usZ3A+SPWKI1s0K7cMXq1YizXc0+yCXgecC1wbFXth94HQJJjlnnOVmArwFOe8pRxlCENxTBXi0beoZrkicBHgd+uqm8P+ryq2lFV81U1Pzc3N2oZkqQ+I4V7ksfRC/YPV9XHutn3Jzmue/w44MBoJUrTZzev9W6Uo2UCfBDYW1Xv6ntoF7Clu78FuGr48iRJwxilc38h8FrgxUlu7P69BNgOnJ3kDuDsbrppi0+a0aHhsX7W7pjVWht6h2pV/QuQZR4+a9jXlSSNbixHy0jrjV21Wme4q0lLhfdD8+7eft5Ir7vU8/vXN8rrS+PitWUkqUGGuzRmyw35uJNV02S4S1KDDHdpBYu7bTtwrQeGuyQ1yKNl1ITVdNLj7Lrt4DWr7NwlqUGGuzSAWRlnn1YNs/B/1WgclhmBvwAaxDBDRuvtRKjlTu7S2rFzl6QGGe7SOjMrQ0SabYa7JDXIMfcBrdexUA1v2jsv+99b4xzDXum9O8p7228Qs8vOXZIaZOe+SnYqmpSV/rLTUvNW6rZ9vx667NylRjy0o/Wxrkq5FvprWmln8CA7i92hPBjDXZIa5LDMAOwSNG6DvKdGWWaQ7hdG/6tUy73mpH5nFv8x+nHsdO5/nZYOnLBzl6QGGe7SBD3W33JdT1Yaz1+87CDzxzF2vtRrTHqfw3oZ8zfcJalBjrn38eJHat1jdbnLvf/Xuksd9NvCOPcftMBwp52vztKo1vJ9P+jwympea7mdu5P8fw76QTPpnbcOy0hSgw7Jzr3/U3vxp6YduzR5i4eDxv2ao772cs9bfP2fpeYv9zrTHvKdWOee5Nwktye5M8m2Sa1HkvRoE+nckxwG/DVwNrAP+FKSXVV16yTWt/gkBHeKSqObpfH3WXnd1V7/Z7llppFRk+rczwDurKqvV9X3gSuAzRNalyRpkVTV+F80eSVwblW9oZt+LfDzVfXmvmW2Alu7yWcAt4+9kJVtBL65ButdrfVQ53qoEaxz3KxzfIap8alVNbfUA5PaoZol5j3iU6SqdgA7JrT+gSTZU1Xza1nDINZDneuhRrDOcbPO8Rl3jZMaltkHnNg3fQJw34TWJUlaZFLh/iXg5CQnJXk8cAGwa0LrkiQtMpFhmao6mOTNwD8ChwGXVtUtk1jXiNZ0WGgV1kOd66FGsM5xs87xGWuNE9mhKklaW15+QJIaZLhLUoOaD/eVLoOQZHOSm5LcmGRPkl+cxTr7lvu5JA925xJM3QDb88wk3+q2541J/ngW6+yWObOr8ZYk/zztGrsaVtqev9e3LW/ufvZHz1iNP5nkH5J8pduWF06zvlXUeVSSj3e/79clOW2N6rw0yYEkNy/zeJL8Zff/uCnJ6UOtqKqa/UdvZ+6/AU8DHg98BThl0TJP5OF9D88GbpvFOvuW+zzwaeCVs1gncCbwyXXwcz8SuBV4Sjd9zCzWuWj5lwGfn7UagT8A/qy7Pwc8ADx+But8J3BJd/+ZwO5p/8y7db8IOB24eZnHXwJ8ht75Qs8Hrh1mPa137iteBqGq/re6LQocwaKTraZk0Ms1vAX4KHBgmsX1WS+XlRikzl8HPlZV9wBU1Vps09Vuz1cDl0+lsocNUmMBP5Ek9JqlB4CD0y1zoDpPAXYDVNVtwKYkx063TKiqa+hto+VsBj5UPV8Ejkxy3GrX03q4Hw/c2ze9r5v3CElekeQ24FPA66dUW78V60xyPPAK4P1TrGuxgbYn8ILuK/pnkpw6ndIeYZA6nw4cleQLSa5P8rqpVfewQbcnSQ4HzqX34T5Ng9T4HuBZ9E5U/Crw1qr64XTK+5FB6vwK8MsASc4AnkrvBMtZM/D74rG0Hu4rXgYBoKo+XlXPBM4H3j7popYwSJ3vBt5WVQ9OvpxlDVLnl+ld7+I5wF8Bn5h0UUsYpM4NwM8C5wHnAH+U5OmTLmyRgd6fnZcB/1pVj9XxTcIgNZ4D3Ag8GXgu8J4kT5psWY8ySJ3b6X2g30jvW/ANTP8bxiBW875YVut/rGNVl0GoqmuS/HSSjVU1zYsMDVLnPHBF75svG4GXJDlYVZ+YSoU9K9ZZVd/uu//pJO+d0e25D/hmVX0H+E6Sa4DnAF+bTok/qmHQ9+cFTH9IBgar8UJgeze8eWeSu+iNaV83nRKBwd+bF0JvpyVwV/dv1ozn8i1rsUNhijsuNgBfB07i4Z0spy5a5md4eIfq6cA3HpqepToXLX8Za7NDdZDt+VN92/MM4J5Z3J70hhF2d8seDtwMnDZrdXbL/SS9MdojZvRn/j7gT7r7x3a/QxtnsM4j6Xb0Ar9Bb1x7qtuzr5ZNLL9D9TweuUP1umHW0XTnXstcBiHJG7vH3w/8CvC6JD8A/g/4teq28IzVueYGrPOVwG8mOUhve14wi9uzqvYm+SxwE/BD4ANVteShaWtZZ7foK4DPVe9bxlQNWOPbgcuSfJVeIL2tpvtNbdA6nwV8KMmD9I6UumiaNT4kyeX0jirbmGQfcAnwuL46P03viJk7ge/SfdtY9Xqm/HsnSZqC1neoStIhyXCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDfp/nSelzv2pY4wAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt \n",
    "plt.hist(stats.norm(0,1).cdf(trace.posterior[\"l_z\"])[0,:,0],bins=200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3.87217648, 4.11651589, 4.16778903, 4.12829874, 3.78085534,\n",
       "        3.82195114, 3.85739836, 4.05119681, 3.68041513, 4.15573983]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lambda_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import arviz as az\n",
    "summarize = az.summary(trace,\"zeta\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.8444946583966609, 3.7516082409319318e-137)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy.stats as stats\n",
    "stats.pearsonr(summarize[\"mean\"].values.flatten(),zeta_p.flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "data = []\n",
    "import time\n",
    "\n",
    "while True:        \n",
    "    try:\n",
    "        try:\n",
    "            with open(\"./indicator.pkl\",\"rb\") as f:\n",
    "                indicator_saved = pickle.load(f)\n",
    "                \n",
    "        except Exception as e:\n",
    "            print(1)\n",
    "            with open(\"./indicator.pkl\",\"wb\") as f:\n",
    "                pickle.dump([],f) \n",
    "                \n",
    "            \n",
    "        with open(\"./indicator.pkl\",\"wb\") as f:\n",
    "            indicator_saved.append(data)\n",
    "            pickle.dump(indicator_saved,f) \n",
    "        break\n",
    "        \n",
    "    except Exception as ie:\n",
    "        print(ie)\n",
    "        time.sleep(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"./indicator.pkl\",\"wb\") as f:\n",
    "        indicator_saved.append(data)\n",
    "        pickle.dump(indicator_saved,f) "
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "d099968071987d153bdd47ebd8318547dae2f0a1c0c5e02dbbb956a381b0529a"
  },
  "kernelspec": {
   "display_name": "py_mc_kernel",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
